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cycles. 



qqMARKUS Brede 1 

1 CSIRO Marine and Atmospheric Research, CSIRO Centre for Complex System Science, F C Pye Laboratory 
GPO Box 1666, Clunies Ross Street Canberra ACT 2601, Australia 

on 



(N 



< 

d 

13 



> 



PACS 05.45.Xt - Synchronization; coupled oscillators 

PACS 89 . 75 . Fb - Structures and organization in complex systems 

PACS 89 . 75 . He - Networks and genealogical trees 

Abstract. - In this paper we investigate the influence of directed motifs on the synchronization 
properties of Kuramoto oscillators on directed networks. Building different types of sparse directed 
small world networks similar to the Watts and Strogatz procedure we establish that feed forward 
loops favour synchronization on directed networks. The paper highlights the importance of local 
network characteristics for synchronization. 
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Introduction. — Starting with Watts and Strogatz 
introduction of the small world paradigm [1] synchroniza- 
tion phenomena on complex networks have found much 
attention in the recent literature. A question of pervad- 
ing interest in the field is how features of the network 
architecture relate to characteristics of the dynamics on 
the network, e.g. the onset of synchronization or the 
achievement of full synchronization. So far much effort 
has been devoted to an exploration of the stability of the 
fully synchronized state via the master stability function 
(MSF) approach [2] which has allowed considerable insight 
about general structural requirements of network topolo- 
gies to synchronize (see, e.g., [3-9]). One important find- 
ing is that the variance of the in-degrce distribution of 
the network is a strong determinant of the synchroniz- 
ablity [5,6,8], the more in-degree homogeneous networks 
being easier to synchronize than more heterogeneous ones. 
Other studies have demonstrated that small distances, low 
clustering, disassortative degree mixing and an even dis- 
tribution of loads on nodes tend to enhance the stability of 
the fully synchronized state. However, while allowing gen- 
eral insights about a broad class of systems, the MSF ap- 
proach is limited to systems of symmetrically coupled [10] 
identical oscillators and properties of the fully synchro- 
nized state. 



Studying the transition to synchronization in the Ku- 
ramoto system [11] in more detail, a number of recent 
studies has shown that the synchronization properties can 



depend on the whole range of coupling strengths [12,13]. 
Networks may exhibit an onset of synchronization for very 
small coupling, but the transition to full synchronization 
may only occur for very large coupling. Insights in these 
studies have generally been obtained by a finite size scaling 
analysis of the synchronization transition [13-15]. Other 
recent studies have explored the role of the oscillator place- 
ment for synchronization [16-18]. 

Notably, most studies of synchronization on networks 
have focussed on undirected networks. In many cases, 
however, real-world networks are inherently directed. 
Thus, they are characterized by a number of features, as, 
e.g., substantially different in- and out-degree sequences, 
a variety of in- and out-degree correlations, or a more sub- 
tle component organization, that are not found in undi- 
rected networks. One might expect that these characteris- 
tics will contribute to the richness of dynamical phenoma 
on such networks. Indeed, this has been confirmed in a 
couple of recent studies [16, 19] that started to explore the 
relationship between structural features of directed net- 
works and synchronization properties. One finding of this 
work is that the component structure of directed networks 
strongly influences synchronization properties [19]. Fur- 
ther, in [16] a connection between various degree correla- 
tions and synchronization has been made. 

Another observation that adds to the importance of 
studying synchronization on directed networks is that the 
synchronized state can generally be achieved for less cou- 
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pling, i.e. less 'cost'. This may be relevant for technologi- 
cal applications. 

In this paper, we build on these insights and explore the 
influence of small subgraphs or motifs (for a definition of 
network motifs see [21]) for a network's synchronization 
properties. Specifically, we study the transition towards 
synchronization in the Kuramoto model on directed small 
world networks that have been constructed in such a way 
that different (directed) motifs are over- or underrepre- 
scnted, while other relevant network characteristics have 
been kept fixed. This allows to make a strong point that 
directed motifs, i.e. features of the local organization of 
directed networks, can strongly influence the transition 
towards synchronization. 

Model Description. — More specifically, we study 
the Kuramoto system 

(f>i = ojj + a djj sm(<j>j - 4>j), (1) 

3 

where the 4>i,i = 0,...,N— 1 are the phases and the u>i 
the native frequencies of the limit-cycle oscillators, a is 
the coupling strength, and the binary coupling matrix 
ay <G {0, 1} represents the interaction network. For all 
the following experiments, native frequencies are always 
drawn uniformly at random from the interval [—1, 1]. Be- 
cause of its simplicity the Kuramoto model is a suitable 
choice for the dynamics. 

Following a very similar procedure to that of Ref. [1], the 
interaction networks are constructed to be small worlds. 
For this, we start with two different directed ring graph 
substrates, depicted in Fig. 1. The ring graph model (a) is 
built from triangles that arc feed forward loops, while the 
model (b) builds a ring graph from cycles. More specifi- 
cally, to construct ring graphs of range r we set 

r 

Ojj = ]P 8j,i-k (2) 
fc=i 

for model (a) and 

r-l 

a Ti = <S,-,i-i + S i-i+k+i (3) 
fe=i 

for model (b) and i,j = 1,...,N = 2n. In Eq.'s (2) and 
(3) 6i j = 1 if i = j and 6ij = if i ^= j denotes the usual 
Kronecker delta and all indices are understood modulo N. 

It is important to note that both graphs are strongly 
connected, have the same in- and out-degree sequences 
and the same correlations between degrees as well as the 
same triangle densities, i.e. the same (undirected) cluster- 
ing coefficients. Note, however, that the average distances 
between nodes are different on both substrate graphs. Es- 
sentially, cycles allow a walker to also step 'backwards', 
leading to substantially smaller average and maximum dis- 
tances cimax ~ N/3 between nodes for the model (b) in 
comparison to the model (a) for which (i^ax ~ N/2. 




Fig. 1: Illustration of the directed substrate ring graphs to con- 
struct directed small worlds, (a) A ring graph of range 2 built 
from feed forward loops and (b) a ring graph of range 2 built 
from cycles. Note, that both graphs are strongly connected 
and have the same in- and out-degree sequences. Small worlds 
built from cycles are smaller than small worlds build from the 
feed forward loop. 

The construction idea of the ring graphs from loops of 
length three can easily be generalized to loops of arbitrary 
length. For looplength I one may define 

r-l 

aji = Sj t i-i + ^2 8j,i-k(i-i) (4) 
fc=i 

for a ring graph built from feed forward loops of length Z, 
i.e. model (a), and 

r-l 

aji = 6j,i-i + ^2 8j,i+k{l-i), (5) 
fe=i 

i,j= 1, N = (Z — l)(r — l)n, for a ring graph built from 
cycles of length /, see model (b). This construction results 
in substrate graphs with loops no shorter than length I. 
One notes that in this more general situation distances de- 
crease with increasing loop length Z, but are still larger for 
model (a) where d^L ~ N/((r - - 1)) m comparison 
to those for model (b), where d£L ~ N/((r-l){l-l) + l), 
Below we restrict ourselves to graphs with range r = 2. 
We have, however, also explored other values of the range 
only to find no qualitiative difference to the general results 
presented. 

Another important point to note is that following [22] 
one can characterize the synchronizablity of local motifs 
by the coupling strength a* for which the average prob- 
ability that the oscillators on the nodes arc synchronized 
exceeds 1/2. Numerical simulations show that compared 
to feed forward loops cycles are the more synchronizable 
motif in this regard. Hence, based on shorter distances 
and higher synchronizability of the local motifs one should 
expect small worlds built from the configuration b) to be 
more synchronizable than those built from configuration 
a). As we will see below, this is however not the case. 
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From the substrate graphs small worlds are obtained in 
the following way. Iteratively, a node and one of its incom- 
ing links arc picked randomly. The link is then rewired to a 
new randomly selected origin node, provided the resulting 
graph is still strongly connected [20] . If the rewiring dis- 
connects the strong component, a different link is picked 
instead. The procedure is repeated pL times, where L de- 
notes the number of links in the graph. For small p there 
are roughly pL shortcut links in the system. As in [1] 
the parameter p is the shortcut density. Increasing p from 
p = gradually interpolates between a regular and a ran- 
dom graph topology. Note, however, that the graph is not 
fully random even for large p since we have made sure of 
the strong connectedness and an equal number of incom- 
ing links for all nodes during the rewiring procedure [23] . 

Results. — The degree of synchronization between the 
oscillators of (1) is usually measured by the order param- 
eter 



r(t) 



l/iV£exp(i^(f)) 



(6) 



If the phases of the N oscillators are desynchronized one 
has r ~ l/ViV while r ~ 0(1) when a macroscopic part 
of the oscillators is synchronized. Increasing the coupling 
strength a between the oscillators one finds a second order 
phase transition. 

In our numerical experiments the order parameter (6) 
is determined by averaging over time (after allowing for 
a sufficient relaxtion time T re i = 100... 400 depending on 
network size), and the randomness in the network, ini- 
tial conditions (j>i(t = 0) which are randomly drawn from 
(— 7T, 7r], and the oscillator ensemble, i.e. 



2T rol 

r=(l/T re] J2 K*)>, 

f=T re] 



(7) 



where brackets indicate the average over the different 
sources of randomness and T Te \ is the relaxation time. 

Figure 2 illustrates some first numerical results for the 
synchronization transition on small world networks of size 
N = 1600 with shortcut density p = 0.05 built from the 
two different motifs. The data have been obtained for 
oscillators chosen uniformly at random from [—1,1] as av- 
erages from a numerical integration of the system (1) over 
200 different randomly constructed small world networks 
built from the substrate graphs introduced above. Clearly, 
the onset of synchronization as well the transition to full 
synchronization are achieved for lower coupling for the 
small worlds built from feed forward cycles. 

To understand the difference in the synchronization 
behaviours between both systems, it is insightful to in- 
vestigate the evolution of the sizes and the structures 
of synchronized components as the coupling strength a 
is increased. Synchronized components can be defined 
based on differences in the average frequencies ZJ = 
limT-+oo 1/T Jq <fi(t)dt. Since the frequency resolution in 



our numerical simulations is given by l/(2T rc i), two os- 
cillators 1 and 2 can be regarded as mutually entrained 
if \W\ — ZZJ2 1 < l/(2T rc i). The synchronized in- or out- 
components of a given oscillator can then be defined as the 
entrained oscillators that can be reached following in- or 
out-links via entrained oscillators. Likewise, synchronized 
strongly connected components (SSCCs) arc maximal sets 
of entrained oscillators that can reach each other by con- 
nected paths via entrained oscillators. In the following we 
denote the average size of SSCCs by -/V sscc and the average 
relative size by n sscc = N sscc /N. 

Panel (b) of Fig. 2 gives the evolution of the relative 
average size of SSCCs when the coupling a is increased. 
Crucially, as one may already expect from the construc- 
tion of the substrate graphs, one notes that N sscc is very 
small almost up to the transition point for small worlds 
constructed by feed forward loops whereas N sscc grows 
gradually with increasing coupling for small worlds con- 
structed from cycles. Moreover, as the distributions of 
the sizes of SSCCs for small coupling displayed in Fig. 2c 
exemplify, there are many trivial and very small SSCCs in 
the first case, whereas many more large SSCCs of roughly 
comparable sizes are found in the second case. Global 
synchronization in the first case can thus be attained by a 
large SSCCs recruiting small or trivial SSCCs, whereas in 
the second case it arises from the competition and merger 
of SSCCs of comparable sizes. Moreover, as we observed 
before, the feed forward motif per se is less synchronizablc 
than the cycle motif. Hence small synchronized compo- 
nents in the first case can be more easily broken up and 
recruited into larger synchronized clusters than in the sec- 
ond case. Superior synchronization on small worlds con- 
structed by feed forward loops is thus attained, because 
the recruitment of small synchronized components into a 
large synchronized component can be achieved with less 
coupling than on small worlds built from cycles, where 
competition between large synchronized components of 
the same size exists. 

Figure 2, however, gives only a first indication of the 
synchronization properties and we support our first result 
by a more comprehensive set of experiments below. A 
more rigorous quantification of the synchronization tran- 
sition can be obtained from a finite size scaling (FSS) anal- 
ysis. In this, we follow the approach of [14, 15] and assume 
the scaling ansatz 



r(<7, N) = N- a F(N^ v (a - a c )) 



(8) 



close to the critical point. In Eq. (8) F(-) is a universal 
scaling function and a and v are critical exponents. The 
exponent v describes the divergence of the correlation vol- 
ume close to the critical point. From 



r(cr) ~(<j- <7 c y 



(9) 



in the thermodynamic limit one obtains a = (3/v. As 
in [14, 15] the critical point can be obtained by plotting 
rN a as a function of the control parameter a for different 
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Fig. 2: (a) Comparison of the transitions to synchronization on small worlds built from feed forward loops and from cycles, 
(b) Average relative size of the synchronized strongly connected components for both systems, (c) Distribution of the sizes of 
strongly connected synchronized components for a = .8. Parameters are N = 1600, T rc i = 300, range r = 2 and looplength 
/ = 3, density of shortcuts p = 0.05. 



system sizes. From (8) one realizes that the curves for dif- 
ferent N intersect at the critical point for the right choice 
of the critical exponent a. To determine the critical points 
we integrate Eq. (1) for several system sizes comprising up 
to N — 6400 nodes. We then proceed with the FSS anal- 
ysis. In Figure 3 an example of such an analysis is carried 
out for small worlds (r = 2, / = 3) with shortcut density 
p = .05 constructed from feed forward loops according to 
model (a) and cycles according to model (b), respectively. 
One finds that for both small world models for a — .25 
plotting rN a vs. the coupling strength a data points for 
all system sizes intersect at a unique crossing point, the 
critical point a c . Having determined the critical coupling 
<7 C and the exponent a one can then again use the FSS as- 
sumption of (8) to make all data points collapse to a single 
smooth curve (cf. insets of the panels of Fig. 3). This al- 
lows the determination of the exponent v as the exponent 
for which the best data collapse when plotting rN a vs. 
(cr — a c )N lL/ is obtained. The analysis yields v = 3.0(1). 
One notes that the value of v for both directed small world 
models is thus different from the value of v determined for 
undirected small worlds in [14]. 

The more detailed analysis yields rr c = .98(2) for model 
(a) while er c = 1.78(2) for model (b). Thus, as one might 
already anticipate from Fig. 2 the synchronization tran- 
sition on small world networks built from feed forward 
cycles occurs for substantially smaller coupling than for 
small world networks built from cycles. This appears par- 
ticularly surprising, since small worlds built by cycles are 
typically much smaller than small worlds constructed from 
feed forward loops. In- and out degree sequences of the 
graphs and the triangle densities being the same, one con- 
cludes that the difference in the onset of synchronization 
is related to the local organization of the networks, i.e. the 
different motifs that they are constructed from. 

This result proves robust if one varies the shortcut den- 
sity of the small worlds, gradually tuning them towards 
random graphs, cf. Fig. 4. For all shortcut densities we 
experimented with the FSS analysis yields ai < ■ 



Naturally, however, for large p as more and more links are 
rewired the local organisation is destroyed and the critical 
points of the two small world models converge. 

To further strengthen the main assertion of the pa- 
per that feed forward loops facilitate synchronization, we 
also investigated a probabilistic variant of the small world 
model, for which the fraction of feed forward loops and cy- 
cles can be systematically tuned. For this, we consider the 
ring graph with r = 2 and I = 3 and set a# = r5 J - l ,_i+<5j l ,_2 
with probability qg and dji = + <5j.i+i otherwise. As 

implied by the notation the parameter qg gives the ratio of 
feed forward loops and cycles. One should, however, note 
that in this model the in-degree sequence of the graphs 
is no longer homogeneous. Instead, the variance of the 
in-degree sequence is cr? = 2qg (1 — qg), which has a max- 
imum for equal fractions of feed forward loops and cycles, 
i.e. for qg = 1/2. Hence, as the in-degree heterogeneity 
strongly influences a network's synchronization properties, 
one has two competing tendencies when qg is tuned. First, 
increasing qg from zero onwards the in-degree variance of 
the graphs grows strongly. Second, increasing qg also in- 
creases the fraction of feed forward cycles. 

Both competing tendencies can be seen from the data in 
Fig. 5, which gives the dependence of the critical coupling 
a c on qg. By the above discussed increase in in-degree 
variance, initially, replacing cycles by feed forward loops 
deteriorates the networks' synchronizabilities. The nega- 
tive impact of the increased heterogeneity competes with 
the positve impact of the enhanced feed forward loop den- 
sity. The data in Fig. 5 verify that o~ c reaches a maximum 
before the in-degree variance reaches its maximum. 

Another point to note is that the in-degree variance is 
the same for graphs built with a given density q = qg (for 
feed forward loops) or q = 1 — qg (for cycles) of one of the 
motifs. Thus, comparing the synchronization threshold 
for graphs built with qg = q and qg = 1 — q excludes the 
effect of different in-degree heterogeneities. For this com- 
parison Fig. 5 clearly demonstrates that larger densities 
of feed forward loops always entail a lower synchronization 
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Fig. 4: Dependence of the critical point of the synchronization 
transition a c on the shortcut density p for small worlds with 
r = 2 and I — 3. The data are obtained from a FSS analysis, 
error bars are smaller than the symbols. Lines are just guides 
for the eye. 



previous results for triangles one again observes that the 
graphs constructed from feed forward loops are more syn- 
chronizablc than the graphs built from cycles. 

The data of Fig. 6 indicate that the critical points for 
both models converge when loops become very long. One 
may surmise that for large I the impact of shorter loops 
introduced by the shortcut links becomes dominant, thus 
obliterating the difference between the models. 



Fig. 3: FSS analysis of the transition to synchronization for 
small world networks with short cut density p — 0.05. (a) The 
transition for small worlds built from feed forward loops, ot — 
.25, <T C = .98(2). (b) The transition for small world networks 
built from cycles, a = .25, a c — 1.78(2). The insets of the 
panels (a) and (b) show the collapse of data for different system 
sizes when plotting Nr a vs. (a — cr c )N 1 ^ 1 ' for v = 3.0(1). Data 
points represent averages over at least 100 independent runs, 
error bars are smaller than the size of the data points. 

threshold a c . 

A natural extension of the conventional small world 
model with immediate nearest neighbour connections is to 
construct the ring graphs from loops of arbitrary length I, 
cf. Eq.'s (4) and (5). In this way the influence of general 
feed forward loops and cycles on the onset of synchro- 
nization can be explored. In Fig. 6 simulation results 
for the dependence of the critical point of small worlds 
with shortcut density p = 0.05 on the loop length are 
displayed. Generally, introducing longer loops reduces av- 
erage distances and thus favours synchronization. This is 
confirmed by a gradual decline of the critical points with 
increasing loop length in Fig. 6. In agreement with the 



Conclusions. — In this paper we have studied the 
synchronization transition in several models of directed 
small world graph topologies. We stress the point that the 
synchronizablity of a small world strongly depends on the 
particulars of the local organization. Generally, having the 
same degree sequences and even though being larger in av- 
erage distances, small world graphs built purely from feed 
forward loops are found to exhibit a transition towards full 
synchronization at a lower critical coupling strength than 
small worlds built purely from cycles. This point is cor- 
roborated by the study of the synchronization threshold of 
small worlds constructed with a tuneable ratio of feed for- 
ward loops and cycles. Other parameters being the same, 
small worlds with larger feed forward loop densities always 
exhibit a lower threshold to synchronization. 

Similar to the spirit of Rcf. [22] we also characterised 
the synchronizability of the standalone motifs. The find- 
ing that the more synchronizable motif gives rise to infe- 
rior synchronization properties globally cautions to con- 
clude about global synchronization from local subgraph 
properties. It thus appears an appealing agenda for fur- 
ther research to study the influence of directed motifs on 
a network's synchronization properties. 
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Fig. 5: Dependence of the critical point of the synchronization Fig. 6: Dependence of the critical point a c on the length of 
transition a c on the ratio of feed forward loops to cycles qg for loops for feed forward loops and cycles for small worlds with 



small worlds with shortcut density p 
analysis, lines are guides for the eye. 



.2. Results from FSS shortcut density p = .05. Results obtained from a FSS analysis. 
Error bars are below the size of the symbols. 
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